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Abstract. A code for the numerical evaluation of hyperelliptic theta-func- 
tions is presented. Characteristic quantities of the underlying Riemann surface 
such as its periods are determined with the help of spectral methods. The 
code is optimized for solutions of the Ernst equation where the branch points 
of the Riemann surface are parameterized by the physical coordinates. An 
exploration of the whole parameter space of the solution is thus only possible 
with an efficient code. The use of spectral approximations allows for an efficient 
calculation of all quantities in the solution with high precision. The case of 
almost degenerate Riemann surfaces is addressed. Tests of the numerics using 
identities for periods on the Riemann surface and integral identities for the 
Ernst potential and its derivatives are performed. It is shown that an accuracy 
of the order of machine precision can be achieved. These accurate solutions are 
used to provide boundary conditions for a code which solves the axisymmetric 
stationary Einstein equations. The resulting solution agrees with the theta- 
functional solution to very high precision. 



1. Introduction 

Solutions to integrable differential equations in terms of theta-functions were in- 
troduced with the works of Novikov, Dubrovin, Matveev, Its, Krichever, . . . (see [7, 
19, 28, 2]) for the Korteweg-de Vries (KdV) equation. Such solutions to e.g. the 
KdV, the Sine-Gordon, and the Non-linear Schrodinger equation describe periodic 
or quasi-periodic solutions, see [8, 2]. They are given explicitly in terms of Riemann 
theta-functions defined on some Riemann surface. Though all quantities entering 
the solution are in general given in explicit form via integrals on the Riemann 
surface, the work with theta-functional solutions admittedly has not reached the 
importance of soliton solutions. 

The main reason for the more widespread use of solitons is that they are given in 
terms of algebraic or exponential functions. On the other hand the parameterization 
of theta-functions by the underlying Riemann surface is very implicit. The main 
parameters, typically the branch points of the Riemann surface, enter the solutions 
as parameters in integrals on the Riemann surface. A full understanding of the 
functional dependence on these parameters seems to be only possible numerically. 
In recent years algorithms have been developed to establish such relations for rather 
general Riemann surfaces as in [33] or via Schottky uniformization (see [2]), which 
have been incorporated successively in numerical and symbolic codes, see [32, 18, 
14, 5, 6] and references therein (the last two references are distributed along with 
Maple 6, respectively Maple 8, and as a Java implementation at [36]). For an 
approach to express periods of hyperelliptic Riemann surfaces via theta constants 
see [9]. 
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These codes axe convenient to study theta- functional solutions of equations of KdV- 
type where the considered Riemann surfaces are 'static', i.e., independent of the 
physical coordinates. In these cases the characteristic quantities of the Riemann 
surface have to be calculated once, just the comparatively fast summation in the 
approximation of the theta series via a finite sum as e.g. in [6] has to be carried out 
in dependence of the space-time coordinates. 

The purpose of this article is to study numerically thcta-functional solutions of the 
Ernst equation [10] which were given by Korotkin [25] . In this case the branch points 
of the underlying hyperelliptic Riemann surface are parameterized by the physical 
coordinates, the spectral curve of the Ernst equation is in this sense 'dynamical'. 
The solutions are thus not studied on a single Riemann surface but on a whole 
family of surfaces. This implies that the time-consuming calculation of the periods 
of the Riemann surface has to be carried out for each point in the space-time. 
This includes limiting cases where the surface is almost degenerate. In addition 
the theta-functional solutions should be calculated to high precision in order to be 
able to test numerical solutions for rapidly rotating neutron stars such as provided 
e.g. by the spectral code LORENE [35]. This requires a very efficient code of high 
precision. 

We present here a numerical code for hyperelliptic surfaces where the integrals 
entering the solution arc calculated by expanding the integrands with a Fast Cosine 
Transformation in MATLAB. The precision of the numerical evaluation is tested 
by checking identities for periods on Riemann surfaces and by comparison with 
exact solutions. The code is in principle able to deal with general (non-singular) 
hyperelliptic surfaces, but is optimized for a genus 2 solution to the Ernst equation 
which was constructed in [22, 23] . We show that an accuracy of the order of machine 
precision 10^^^) can be achieved at a space-time point in general position with 32 
polynomials and in the case of almost degenerate surfaces which occurs e.g., when 
the point approaches the symmetry axis with at most 256 polynomials. Global 
tests of the numerical acciiracy of the solutions to the Ernst equation are provided 
by integral identities for the Ernst potential and its derivatives: the cqiiality of 
the Arnowitt-Deser-Misner (ADM) mass and the Komar mass (see [24, 34]) and 
a generalization of the Newtonian virial theorem as derived in [15]. We use the 
so determined numerical data for the theta-functions to provide 'exact' boundary 
values on a sphere for the program library LORENE [35] which was developed for a 
numerical treatment of rapidly rotating neutron stars. LORENE solves the boundary 
value problem for the stationary axisymmetric Einstein equations with spectral 
methods. We show that the theta-functional solution is reproduced to the order of 
10"" and better. 

The paper is organized as follows: in section 2 we collect useful facts on the Ernst 
equation and hyperelliptic Riemann surfaces, in section 3 we summarize basic fea- 
tures of spectral methods and explain our implementation of various quantities. 
The calculation of the periods of the hyperelliptic surface and the non-Abelian line 
integrals entering the solution is performed together with tests of the precision of 
the numerics. In section 4 we check integral identities for the Ernst potential. The 
test of the spectral code LORENE is presented in section 5. In section 6 we add some 
concluding remarks. 

2. Ernst equation and hyperelliptic Riemann surfaces 

The Ernst equation for the complex valued potential 8 (we denote the real and the 
imaginary part of £ with / and h respectively) depending on the two coordinates 
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(p, C) can be written in the form 



1 



(1) [Spp + -Sp + Sccj =^p+ ^i- 

The equation has a physical interpretation as the stationary axisymmetric Einstein 
equations in vacuum (sec appendix and references given therein). Its complete 
integrability was shown by Maison [29] and Belinski-Zakharov [1]. For real Ernst 
potential, the Ernst equation reduces to the axisymmetric Laplace equation for 
\n£. The corresponding solutions are static and belong to the so called Weyl class, 
see [27]. 

Algebro-geometric solutions to the Ernst equation were given by Korotkin [25] . The 
solutions are defined on a family of hyperelliptic surfaces jC{^,^) with ^ = ( — ip 
corresponding to the plane algebraic curve 

(2) „' = {K-0{K-of[iK-E,){K-Fi), 

i=l 

where g is the genus of the surface and where the branch points Ei, Fi are inde- 
pendent of the physical coordinates and for each n subject to the reality condition 

E„ = FnOr En,Fr,GR. 

Hyperelliptic Riemann surfaces are important since they show up in the context of 
algebro-geometric solutions of various integrable equations as KdV, Sine-Gordon 
and Ernst. Whereas it is a non-trivial problem to find a basis for the holomorphic 
differentials on general surfaces (see e.g. [5]), it is given in the hyperelliptic case 
(see e.g. [2]) by 

fdK KdK Ks-^dK 

(3) di>k={ , 

VMM M 

which is the main simplific;ation in the use of these surfaces. We introduce on £ a 
canonical basis of cycles (ak,bk), k = 1, . . . ,n. The holomorphic differentials dwfc 
are normalized by the condition on the a-periods 



(4) / 

J a 



du>k = 2TriSu 



The matrix of 6-pcriods is given by B^fe = dojk- The matrix B is a so-called 
Riemann matrix, i.e. it is symmetric and has a negative definite real part. The 
Abel map a; : £ — > Jac(£) with base point E\ is defined as a;(P) = duj^, where 
Jac(£)is the Jacobian of C The theta-function with characteristics corresponding 
to the curve C is given by 

(5) epq(x|B)= ^ exp|i(B(p + n),(p + n)) + (p + n,2i7rq + x)|, 

where x e is the argument and p, q G are the characteristics. We will 
only consider half-integer characteristics in the following. The theta-function with 
characteristics is, up to an exponential factor, equivalent to the theta-function with 
zero characteristic (the Riemann theta-function is denoted with G) and shifted 
argument, 

(6) epq(x|B) = e(x Bp 2i7rq) exp j i(Bp, p) + (p, 2i7rq x) I . 



We denote by diopq a differential of the third kind, i.e., a 1-form which has poles in 
P,Q € C with respective, residues and —1. This singularity structure charac- 
terizes the differentials only up to an arbitrary linear combination of holomorphic 
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differentials. The meromorphic differentials can be normalized by the condition 
that all a-periods vanish. We use the notation oo^ for the infinite points on dif- 
ferent sheets of the curve £, namely ^/K^^^ ±1 as if ^ oo^. The differential 
<^cx3+oo- is given up to holomorphic differentials by —K^dK/fj,. It is well known 
that the 6-periods of normalized differentials of the third kind can be expressed in 
terms of the Abel map (see e.g. [8]), 

(7) / dojpQ = uJk{P) - Wfc((5), k = l,...,g. 



In [20, 21] a physically interesting subclass of Korotkin's solution was identified 
which can be written in the form 

.ox c_ epqi-.-(^ + )+u) , 

^ ' epq(u;(oo-) + u) ' 

where u = {uk) G and where 

(9) I=^.J^\nG{K)du;^+^-{K), uu = ^. j^\nG{K) duju- 

r is a piece-wise smooth contour on C and G{K) is a non-zero Holder-continuous 

function on T . The contour V and the function G have to satisfy the reality condi- 
tions that with if e r also K and G{K) = G{K); both are independent of the 
physical coordinates. 

In the following we will discuss the example of the solution constructed in [22, 23] 

which can be interpreted as a disk of coUisionless matter. For a physical interpre- 
tation sec [13]. The solution is given on a surface of the form (2) with genus 2. 
The branch points independent of the physical coordinates are related through the 
relations Ei = Fi, i = 1,2 and Ei = —E2. The branch points are parameterized by 
two real parameters A and 6. Writing Ef = a -\- i(i with real a, (3, we have 



5 1 52 

(10) " = -1+2' ^=Va^+'^"T- 

The contour T is the piece of the covering of the imaginary ax;is in the upper sheet 
between [— i,i], the function G has the form 



(11) GiK) = V^^^±^. 

The physical parameters vary between ^ = 0, the solution which was first given 
in [30], and 5,, = 2(1 + y^l + l/A^), the static Hmit in which /? = 0. In the 
latter case the Riemann surface degenerates, the resulting Ernst potential (8) is 
real and be expressed in terms of objects corresponding to the surface Co of genus 
defined by the relation /Iq = {K — £,){K — The parameter A varies between 
A = 0, the so-called Newtonian limit where the branch points Ei, Fi tend to infinity. 
Since G is also of order A in this limit, the lowest order contributions are again 
real and defined on the surface Cq. This case corresponds to the disk limit of 
the Maclaurin ellipsoids, see [3]. The upper limit for A is infinity for (5 ^ and 
Ac = 4.629 ... for ^ = 0. The limiting situation is special in the second case since the 
resulting spacetime is no longer asymptotically flat and since the axis is singular. 
The invariant circumference of the disk is zero in this case which implies that the 
disk shrinks to a point for an observer in the exterior of the disk, see [13]. 

For physical reasons the solution was discussed in [13] in dependence of two other 
real parameters e and 7. Here e is related to the redshift of photons emitted at the 
center of the disk and detected at infinity. It varies between in the Newtonian 
limit, and 1 in the ultra-relativistic limit, where photons cannot escape to infinity. 
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Thus, e is a measure of how relativistic the situation is. The parameter 7 is a 
measure of how static the solution is, it varies between 0, indicating the static hmit 
and 1. For the functional relations between e, 7 and A, S see [13]. The constant 
(with respect to the physical coordinates) to appear in the following can be 
considered as a natural scale for the angular velocities in the disk, for a definition 
see [13]. 

The coordinate p can take all non- negative real values, the coordinate C all real 
values. The example we are studying here has an equatorial symmetry, 

(12) £ip,-0=£{p,0. 

It is therefore sufficient to consider only non-negative values of The case p = 
corresponds to the axis of symmetry where the branch cut ^] degenerates to a 
point. As was shown in [21, 13], the Ernst potential can be written in this limit in 
terms of theta-fimctions on the elliptic surface £1 defined by pf — [K"^ — of' -\- 0^ , 
i.e. the surface C with the cut [5, f] removed. Near the axis the Ernst potential has 
the form (see [11, 21]) 

(13) £:(p,C)=^o(C)+p'^i(C)+0(p'); 

here £q and £\ are independent of p, £0 is the axis potential. This formula could 

be used to calculate the potential close to the axis. However we considered only 
values of p greater than 10~^ and did not experience any numerical problems. 
Consequently we did not use formula (13). 

For large values of r = |^], the Ernst potential has the asymptotic expansion 

(H) . = l-|!,?^_?!|<,0,l/.3); 

here the constants (with respect to ^) m and J arc the ADM-mass and, respectively, 
the angular momentum of the space-time. They can be calculated on the axis in 
terms of elliptic theta-functions, see [13]. Formula (14) is used for values of r > 10^. 

In the limit ^ = E2, the Ernst potential can be given on the surface Sq of genus 
obtained by removing the cuts [^, £] and [E2, F2] from the surface C. The potential 
can thus be given in this case in terms of elementary functions, see [13]. 

In the equatorial plane C = 0, the Riemann surface C has an additional involution 
K — > —K as can be seen from (2). This implies that the surface can be considered 
as a covering of an elliptic surface, see [2, 21]. The theta-functions in (8) can be 
written as sums of theta-functions on the covered surface and on the Prym variety 
which happens to be an elliptic surface as well in this case. We use this fact at the 
disk (C = 0, p < 1), where the moving branch points are situated on F. There all 
quantities can be expressed in terms of quantities defined on the Prym surface T,yj 
defined by nl = {K + p^){{K - af + f3'^), see [13]. 

3. Numerical implementations 

The numerical task in this work is to approximate and evaluate analytically defined 

functions as accurately and efficiently as possible. To this end it is advantageous to 
use (pseudo-)spectral methods which are distinguished by their excellent approxi- 
mation properties when applied to smooth functions. Here the functions are known 
to be analytic except for isolated points. In this section we explain the basic ideas 
behind the use of spectral methods and describe in detail how the theta-functions 
and the Ernst potential can be obtained to a high degree of accuracy. 
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3.1. Spectral approximation. The basic idea of spectral methods is to approxi- 
mate a given function / globally on its domain of definition by a linear combination 

N 

f ~ y^Qfc^fc, 

fe=0 

where the function are taken from some class of functions which is chosen ap- 
propriately for the problem at hand. 

The coefficients (ik arc determined by requiring that the linear combination should 
be 'close' to /. Thus, one could require that \\.f — '^^^q afe<^fc|| should be minimal for 
some norm. Another possibility is to require that ^/ — X^^g ^k4>k^ Xi^ = ' = 
: A'' with an appropriate inner product and associated orthonormal basis xi- This 
is called the Galcrkin method. Finally, one can demand that f{xi) = X]^o '^■k<Pk{xi) 
at selected points {xi)i^o:N- This is the so called collocation method which is the 
one we will use in this paper. In this case the function values fi = f{xi) and the 
coefficients ak are related by the matrix = (j)k{xi). 

The choice of the expansion basis depends to a large extent on the specific problem. 
For periodic functions there is the obvious choice of trigonometric polynomials 
(t>k{x) = ex.p{2-jTik/N) while for functions defined on a finite interval the most 
used functions arc orthogonal polynomials, in particular Chcbyshcv and Legcndre 
polynomials. While the latter are important because of their relationship with the 
spherical harmonics on the sphere, the former are used because they have very good 
approximation properties and because one can use fast transform methods when 
computing the expansion coefficients from the function values provided one chooses 
the collocation points xi = cos{'itI/N) (see [12] and references therein). We will use 
here collocation with Chebyshev polynomials. 

Let us briefly summarize here their basic properties. The Chebyshev polynomials 
Tn{x) are defined on the interval / = [—1, 1] by the relation 

T„(cos(t)) = cos(nt), where x = cos{t), t G [0,7r]. 

They satisfy the differential equation 

(15) (1 - x^) <l>"{x) - x<t>'(x) + n^(t>(x) = 0. 

The addition theorems for sine and cosine imply the recursion relations 

(16) T„+i(x) - 2xTn{x)+Tn-i{x) = 0, 
for the polynomials T„ and 

(17) I^l±M-IkzM=2T^{x) 

for their derivatives. The Chebyshev polynomials are orthogonal on I with respect 
to the hermitian inner product 

{f,9) = f_J{x)-9{x) '^'^ 



We have 

TT 

(-18) ^mi^n) ~ Cm~^ ^mn 

where cq = 2 and Q = 1 otherwise. 

Now suppose that a function / on 7 is sampled at the points xi = cos^nl/N) and 
that J2n=o ^nTn is the interpolating polynomial. Defining cq = Cjv = 2, c„ = 1 for 
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< n < A'' in the discrete case and the numbers F„ = c„a„ we have 

N N 
fl = Y^ (^nTnixi) = J2 anTn{cOs{TTl/N)) 



n=0 n=0 
N N 



= a„ cos{Trnl/N) = — cos{'Knl/N) 

This looks very much hke a discrete cosine series and in fact one can show [4] that 
the coefBcients F„ are related to the values /; of the function by an inverse discrete 
Fourier transform (DCT) 

2 ^ f 
Fn = ^,y2-cos{Trnl/N). 

iV ^—^ Ci 

1=0 ' 

Note, that up to a numerical factor the DCT is idcmpotcnt, i.e., it is its own inverse. 
This relationship between the Chebyshev polynomials and the DCT is the basis for 
the efficient computations because the DCT can be performed numerically by using 
the fast Fourier transform (FFT) and pre- and postprocessing of the coefficients [12]. 
The fast transform allows us to switch easily between the representations of the 
function in terms of its sampled values and in terms of the expansion coefBcients 
an (or F„). 

The fact that / is approximated globally by a finite sum of polynomials allows us to 
express any operation applied to / approximately in terms of the coefficients. Let us 
illustrate this in the case of integration. So we assume that f = Pn = '^n=o ^nTn 
and we want to find an approximation of the integral for p^, i.e., the function 



F{x) = £j{s)ds 



so that F'{x) = f{x). We make the ansatz F{x) = J2n=o^nTn{x) and obtain the 
equation 

N N 

F' = J2bnK = Y,<^nTn = f. 
n=0 n=0 

Expressing r„ in terms of the using (17) and comparing coefiicients implies the 
equations 

2ao — a2 a„_i — a„+i ^ „ ^ ^ ,r i. ojv-i 
b, = ^—, bn = forO<n<7V, = 

between the coefficients which determines all bi in terms of the a„ except for 60- 
This free constant is determined by the requirement that F[—l) = which implies 
(because T„(-l) = (-1)") 

N 
n=l 

These coefficients 6„ determine a polynomial of degree N which approximates 
the indefinite integral F{x) of the iV-th degree polynomial /. The exact function is 
a polynomial of degree N+1 whose highest coefficient is proportional to the highest 
coefficient a^r of /. Thus, ignoring this term we make an error whose magnitude is 
of the order of |aAr| so that the approximation will be the better the smaller {gnI is. 
The same is true when a smooth function / is approximated by a polynomial pn- 
Then, again, the indefinite integral will be approximated well by the polynomial 
qn whose coefficients are determined as above provided the highest coefficients in 
the approximating polynomial pjv are small. 
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From the coefficients 6„ we can also find an approximation to the definite integral 
X^i f{s) ds = F{1) by evaluating 

9Jv(1) = = 2 ^ b2i+i. 

n=0 1=0 

Thus, to find an approximation of the integral of a function / we proceed as de- 
scribed above, first computing the coefficients a„ of /, computing the 6„ and then 
calculating the sum of the odd coefficients. 



3.2. Implementation of the square-root. Tlie Riemann surface C is defined by 
an algebraic curve of the form 

l^^ = {K- 0{K - 1[{K - Ei){K - E,), 

where in our case we have g = 2 throughout. In order to compute the periods and 
the theta-functions related to this Riemann surface it is necessary to evaluate the 
square-root iJ?{K) for arbitrary complex numbers K. In order to make this a 
well defined problem we introduce the cut-system as indicated in Fig. 1. On the 



cutsystem.eps 



Figure 1. Canonical cycles (Pq = C)- 



cut surface the square-root n{K) is defined as in [17] as the product of square-roots 
of monomials 

(19) M = ^/K-EiVK-Ei. 

The square-root routines such as the one available in MATLAB usually have their 
branch-cut along the negative real axis. The expression (19) is holomorphic on the 
cut surface so that we cannot simply take the builtin square-root when computing 
y?{K). Instead we need to use the information provided by the cut-system to 
define adapted square-roots. 

Let arg(z) be the argument of a complex number z with values in ] — 7r,7r[ and 
consider two factors in (19) such as 

^K-P^^/K-P2 
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where Pi and P2 arc two branch-points connected by a branch-cut. Let a = 
a.Yg(P2 — Pi) bo the argument at the Une from Pi to P2. Now we define the 
square-root with branch-cut along the ray with argument a by computing for 
each z e C the square-root s := ^/z with the available MATLAB routine and then 
putting 

s a/2 < arg(s) < a/2 + ir 
— s otherwise 
With this square-root we compute the two factors 

It is easy to see that this expression changes sign exactly when the branch-cut 
between Pi and P2 is crossed. We compute tlie expression (19) by multiplying the 

pairs of factors which correspond to tlic branch-cuts. 

This procedure is not possible in the case of the non-linear transformations we 
are using to evaluate the periods in certain limiting cases. In these cases the root 
is chosen in a way that the integrand is a continuous function on the path of 
integration. 

3.3. Numerical treatment of the periods. The quantities entering formula (8) 
for the Ernst potential are the periods of the Riemann surface and the line integrals 
u and /. The value of the theta-function is then approximated by a finite sum. 

The periods of a hyperelliptic Riemann surface can be expressed as integrals be- 
tween branch points. Since we need in our example the periods of the holomorphic 
differentials and the differential of the third kind with poles at 00^, we have to 
consider integrals of the form 

(20) / n = 0,l,2, 
where the Pi, i, j = I, . . . ,6 denote the branch points of C. 

In general position we use a linear transformation of the form K = ct + d to 
transform the integral (20) to the normal form 

(21) f_^^^±^^^Hit)dt, 

where the ai are complex constants and where H{t) is a continuous (in fact, an- 
alytic) complex valued function on the interval [—1,1]. This form of the integral 
suggests to express the powers t" in the numerator in terms of the first three Cheby- 
shev polynomials To(t) = 1, Ti(t) =t and 12 (i) = 2t^ — 1 and to approximate the 
function H{t) by a linear combination of Chebyshev polynomials 

n>0 

The integral is then calculated with the help of the orthogonality relation (18) of 

the Chebyshev polynomials. 

Since the Ernst potential has to be calculated for all p, C G , it is convenient to 
use the cut-system (1). In this system the moving cut does not cross the immovable 
cut. In addition the system is adapted to the symmetries and reality properties of 
jC. Thus the periods 02 and 62 are related to oi and bi via complex conjugation. 
For the analytical calculations of the Ernst potential in the limit of collapsing cuts, 
we have chosen in [21] cut systems adapted to the respective situation. In the limit 
^ ^ ^ we were using for instance a system where 02 is the cycle around the cut 
[^,^]. This has the effect that only the 6-period 62 diverges logarithmically in this 
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case whereas the remaining periods stay finite as p tends to 0. In the cut systems 
1, all periods diverge as In p. Since the divergence is only logarithmical this does 
not pose a problem for values of p > 10~^. In addition the integrals which have 
to be calculated in the evaluation of the periods are the same in both cut-system. 
Thus there is no advantage in using different cut systems for the numerical work. 

To test the numerics we use the fact that the integral of any holomorphic differential 
along a contour surrounding the cut [Ei, Fi] in positive direction is equal to minus 
the sum of all a-periods of this integral. Since this condition is not implemented in 
the code it provides a strong test for the numerics. It can be seen in Fig. 2 that 
16 to 32 polynomials arc sufficient in general position to achieve optimal accuracy. 
Since MATLAB works with 16 digits, machine precision is in general limited to 14 
digits due to rounding errors. These rounding errors are also the reason why the 
accuracy drops slightly when a higher number of polynomials is used. The use of a 
low number of polynomials consequently does not only require less computational 
resources but has the additional benefit of reducing the rounding errors. It is 
therefore worthwhile to reformulate a problem if a high number of polynomials 
would be necessary to obtain optimal accuracy. These situations occur in the 



apererr . eps 



Figure 2 . Test of the numerics for the a-periods at several points 
in the space-time. The error is shown in dependence of the number 
N of Chebychev polynomials. 

calculation of the periods when the moving branch points almost coincide which 
happens on the axis of symmetry in the space-time or at spatial infinity. As can 
be seen from Fig. 2, for p = 10^"^ and C = 10^ not even 2048 polynomials (this 
is the limit due to memory on the low end computers we were using) produce 
sufficient accuracy. The reason for these problems is that the function H in (21) 
behaves like 1/y/t + p near t = 0. For small p this behavior is only satisfactorily 
approximated by a large number of polynomials. We therefore split the integral in 
two integrals between F2 and {F2 + ^)/2 and between {F2 + and The first 
integral is calculated with the Chcbyshev integration routine after the substitution 
t = \JK — F2 ■ This substitution leads to a regular integrand also at the branch 
point F2. The second integral is calculated with the Chebyshev integration routine 
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after the substitution K — C, = psinh(t). This takes care of the almost collapsing 
cut [5,^]. It can be seen in Fig. 2 that 128 polynomials are sufficient to obtain 
machine precision even in almost degenerate situations. 

The cut-system in Fig. 1 is adapted to the limit ^ ^ in what concerns the 
a-periods, since the cut which collapses in this limit is encircled by an a-cycle. 
However there will be similar problems as above in the determination of the b- 
pcriods. For ^ ^ wc split the integrals for the 6-pcriods as above in two integrals 
between Fi and 0, and and ¥2- For the first integral we use the integration variable 
t = \/K — Fi, for the second K = JfF2 — iQF2 sinht. Since the Riemann matrix 
(the matrix of 6-pcriods of the holomorphic differentials after normalization) is 
symmetric, the error in the numerical evaluation of the 6-periods can be estimated 
via the asymmetry of the calculated Riemann matrix. We define the function 
err{p, C) as the maximum of the norm of the difference in the a-pcriods discussed 
above and the difference of the off-diagonal elements of the Riemann matrix. This 
error is presented for a whole space-time in Fig. 3. The values for p and ( vary 
between 10"'* and 10"*. On the axis and at the disk we give the error for the elliptic 
integrals (only the error in the evaluation of the a-periods, since the Riemann 
matrix has just one component). For ^ — > 00 the asymptotic formulas for the Ernst 
potential are used. The calculation is performed with 128 polynomials, and up to 
256 for 1^1 > 10"^. It can be seen that the error is in this case globally below 10"-^^. 



errorgrey . eps 



Figure 3. A measure for the error in the determination of the 
periods in dependence of the physical coordinates. For p, ^ > 1 we 
use l/p, 1/C as coordinates. 



3.4. Numerical treatment of the line integrals. The line integrals u and / 
in (8) are linear combinations of integrals of the form 



In general position, i.e. not close to the disk and A small enough, the integrals 
can be directly calculated after the transformation K = it with the Chebyshev 



(22) 




I = 0,1,2. 
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integration routine. To test the numerics we consider the Newtonian hmit (A 0) 
where the function InG is proportional to 1 + K'^, i.e. we calculate the test integral 



We compare the numerical with the analytical result in Fig. 4. In general position 
machine precision is reached with 32 polynomials. 



lineerr . eps 



Figure 4. Error in the integrals for the Maclaurin solution in 
dependence of the number N of Chebychev polynomials. 

When the moving cut approaches the path F, i.e., when the space-time point comes 

close to the disk, the integrand in (23) develops cusps near the points ^ and ^. In 
this case a satisfactory approximation becomes difficult even with a large number 
of polynomials. Therefore we split the integration path in [—i, —ip], [—ip,ip] and 
[ip, i]. Using the reality properties of the integrands, we only calculate the integrals 
between and ip, and between ip and i. In the first case we use the transformation 
K = ( + p sinht to evaluate the integral with the Chebyshev integration routine, in 
the second case we use the transformation t = \J K — S^. It can be seen in figure 4 
that machine precision can be reached even at the disk with 64 to 128 polynomials. 
The values at the disk are, however, determined in terms of elliptic functions which 
is more efhcient than the hyperelliptic formulae. 

To treat the case where 5}? is not small, it is convenient to rewrite the function G 

in (11) in the form 



In the limit 5)? oo with 5 finite, the second term in (24) becomes singular 
for K = Q. Even for large but finite, the approximation of the integrand by 
Chebyshev polynomials requires a huge number of coefficients as can be seen from 
Fig. 5. It is therefore sensible to 'regularize' the integrand near K = Q. We consider 



(23) 




(24) 
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instead of the function ln{js — 6 K'^)F{K) where F{K) is a C°° function near K = 0, 
the function 

(25) In (^-^ - SK^^ (^F{K) - F{0) - F'{0)K - ... - -^F(")(0)if" 

The parameter n is chosen such that the spectral coefficients of (25) are of the 
order of 10^^^ for a given number of polynomials, see Fig. 5. There we consider 
the integral 

(26) 



/* 

J —I 



\nG{K)dK 
V(^2-a)2+/?2' 



which has to be calculated on the axis. We show the absolute values of the coeffi- 
cients afe in an expansion of the integrand in Chebyshev polynomials, '^^^^i '^kTk- 
It can be seen that one has to include values of n = 6 in (25). The integral 
Jp In G{K)F{K) is then calculated numerically as the integral of the function (25), 
the subtracted terms are integrated analytically. In this way one can ensure that 



logreg.eps 



Figure 5. Spectral coefficients for the integral (26) for 5 = 1 and 
A = 10^^ in dependence of the number of Chebychev polynomials. 



the line integrals are calculated in the whole space-time with machine precision: 
close to the Newtonian limit, we use an analytically known test function to check 
the integration routine, for general situations we check the quality of the approx- 
imation of the integrand by Chebyshev polynomials via the spectral coefficients 
which have to become smaller than 10~^^. 



3.5. Theta-functions. The theta series (5) for the Riemann theta-function (the 
theta function in (5) with zero characteristic, theta functions with characteristic 
follow from (6)) is approximated as the sum 

^ ^ fl 1 1 

(27) e(x|B) = ^ ^ exp j -n\Bii + nin2Bi2 + -B22 + riiXi + 712X2 > . 
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The value of N is determined by the condition that terms in the series (5) for n > N 
are strictly smaller than some threshold value e which is taken to be of the order 
of 10~^^. To this end we determine the eigenvalues of B and demand that 



where Bmax is the real part of the eigenvalue with maximal real part (B is negative 
definite). For a more sophisticated analysis of theta summations see [6]. In general 
position we find values of N between 4 and 8. For very large values of ( close to the 
axis, can become larger that 40 which however did not lead to any computational 
problems. To treat more extreme cases it could be helpful to take care of the fact 
that the eigenvalues of B can differ by more than an order of magnitude in our 
example. In these cases a summation over an ellipse rather than over a sphere in 
the plane (711,7x2), i-e. different limiting values for rii and 712 as in [6] will be more 
efiicient. 

In our case the computation of the integrals entering the theta-functions was how- 
ever always the most time consuming such that an optimization of the summation 
of the theta-function would not have a noticeable effect. Due to the vectorization 
techniques in MATLAB, the theta summation always took less than 10 % of the 
calculation time for a value of the Ernst potential. Between 50 and 70 % of the 
processor time are used for the determination of the periods. On the used low-end 
PCs, the calculation time varied between 0.4 and 1.2s depending on the used num- 
ber of polynomials. We show a plot of the real part of the Ernst potential for A = 10 
and 5 = 1 in Fig. 6. For p, C > 1, we use l/p, l/C^ as coordinates which makes it 
possible to plot the whole space-time in Weyl coordinates. The non-smoothness 
of the coordinates across p = 1 = 1/ p and C, = 1 = 1/Cis noticeable in the plot. 
Asymptotically the potential is equal to 1. The disk is situated in the equatorial 
plane between p = and p = 1. At the disk, the normal derivatives of / are 
discontinuous. 



fc_10_l.eps 



Figure 6. The real part of the Ernst potential for A = 10 and 
5 = 1 in dependence of the physical coordinates. For p, 1^ > 1 we 
use 1/C as coordinates. 



(28) 




*max 
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The imaginary part of the Ernst potential in this case is given in Fig. 7. It vanishes 
at infinity and at the regular part of the equatorial plane. At the disk, the potential 
has a jump. 



bc_10_l.eps 



Figure 7. The imaginary part of the Ernst potential for A = 10 
and 5 = 1 in dependence of the physical coordinates. For p, ^ > 1 
we use 1/C as coordinates. 



4. Integral identities 

In the previous section we have tested the accuracy of the numerics locally, i.e. at 
single points in the space-time. Integral identities have the advantage that they 
provide some sort of global test of the numerical precision since they sum up the 
errors. In addition they require the calculation of the potentials in extended regions 
of the space-time which allows to explore the numerics for rather general values of 
the physical coordinates. 

The identities we are considering in the following are the well known equivalence of 
a mass calculated at the disk (the Komar mass) and the ADM mass determined at 
infinity, sec [24, 34], and a generalization of the Newtonian virial identity, see [15] 
and the appendix. The derivatives of the Ernst potential occurring in the integrands 
can be related to derivatives of theta-functions, see [21]. Since we are interested 
here in the numerical treatment of theta-functions with spectral methods, we de- 
termine the derivatives with spectral methods, too (see section 3). The integrals 
are again calculated with the Chebyshev integration routine. The main problem in 
this context is the singular behavior of the integrands e.g. at the disk which is a 
singularity for the space-time. As before this will lead to problems in the approx- 
imation of these terms via Chebyshev polynomials. This could lead to a drop in 
accuracy which is mainly due to numerical errors in the evaluation of the integrand 
and not of the potentials which we want to test. An important point is therefore 
the use of integration variables which are adapted to the possible singularities. 

4.1. Mass equalities. The equality between the ADM mass and the Komar mass 
provides a test of the numerical treatment of the elliptic theta-functions at the 
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disk by means of the elliptic theta-functions on the axis. Since this equality is not 
implemented in the code, it provides a strong test. 

The Komar mass at the disk is given by formula (42) of the appendix. In the 

example wc arc considering here, the normal derivatives at the disk can be expressed 
via tangential derivatives (see [13]) which makes a calculation of the derivatives 
solely within the disk possible. We implement the Komar mass in the form 

(29) ruK = t dp , f / + ?(p' - a'f)) ■ 

The integrand is known to vanish as \/\ — at the rim of the disk, which is 
the typical behavior for such disk solutions. Since ^J\-~(P- is not analytic in p, 
an expansion of the integrand (29) in Chebyshcv polynomials in p would not be 
efficient. Wc will thus use t — ^\ — as the integration variable. This takes care 
of the behavior at the rim of the disk. Since in general the integrand in 29 depends 
on this variable can be used in the whole disk. In the ultra-relativistic limit 
for 5 7^ 0, the function / vanishes as p. In such cases it is convenient either to 
take two domains of integration or to use; a different variable of integration. We 
chose the second approach with p = sin a; (this corresponds to the disk coordinates 

(30) ). Yet, strongly relativistic situations still lead to problems since / vanishes in 
this case at the center of the disk as does hp which leads to a '0/0' limit. In Fig. 8 
one can see that the masses are in general equal to the order of 10~^^. In these 
calculations 128 up to 256 polynomials were used. We show the dependence for 
7 = 0.7 and several values of f, as well as for e = 0.8 and several values of 7. The 
accuracy drops in the strongly relativistic, almost static situations (e close to 1, 7 
close to zero) since the Riemann surface is almost degenerate in this case (/? ^ 0). 
In the idtra-rclativistic limit for (5 = 0, the situation is no longer asymptotically flat 
which implies that the masses formally diverge. For e = 0.95, the masses are still 
equal to the order of 10~^^. Not surprisingly the accuracy drops for e = 0.9996 to 
the order of 10~*. 

4.2. Virial-type identities. Generalizations of the Newtonian virial theorem are 

used in numerics (see [15]) as a test of the quality of the numerical solution of the 
Einstein equations. Since they involve integrals over the whole space-time, they 
test the numerics globally and thus provide a valid criterion for the entire range of 
the physical coordinates. 

The identity which is checked here is a variant of the one given in [15] which is 
adapted to possible problems at the zeros of the real part of the Ernst potential, the 
so-called ergosphere, see [13] for the disk solutions discussed here. Eq. (45) relates 
integrals of the Ernst potential and its derivatives over the whole space-time to 
corresponding integrals at the disk. Since the numerics at the disk has been tested 
above, this provides a global test of the evaluation of the Ernst potential. As before, 
derivatives and integrals will be calculated via spectral methods. 

The problem one faces when integrating over the whole space-time is the singular 
behavior of the fields on the disk which represents a discontinuity of the Ernst 
potential. The Weyl coordinates in which the solution is given are not optimal 
to describe the geometry near the disk. Hence a huge number of polynomials is 
necessary to approximate the integrands in (45). Even with 512 polynomials for 
each coordinate, the coefficients of an expansion in Chebyshev polynomials did not 
drop below 10~^ in more relativistic situations. Though the computational limits 
are reached, the identity (45) is only satisfied to the order of 10~^ which is clearly 
related to the bad choice of coordinates. 
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Figure 8. The relative difference of the ADM mass and the Ko- 
mar mass for 7 = 0.7 and several values of e, and for e = 0.8 and 
several values of 7. 



We therefore use for this calculation so-called disk coordinates 77, 9 (see [3]) which 
are related to the Weyl coordinates via 

(30) p + iC = cosh{r] + ie). 

The coordinate rj varies between ry = 0, the disk, and infinity, the coordinate 
between — 7r/2 and 7r/2. The axis is given by ±7r/2, the equatorial plane in the 
exterior of the disk by ^ = and r] ^ 0. Because of the equatorial symmetry, we 
consider only positive values of 6. The surfaces of constant r) arc confocal ellipsoids 
which approach the disk for small 77. For large rj, the coordinates are close to 
spherical coordinates. 

To evaluate the integrals in (45), we perform the ?7-integration up to a value rjo 
as well as the ^-integration with the Chebyshev integration routine. The param- 
eter 770 is chosen in a way that the deviation from spherical coordinates becomes 
negligible, typically r]o = 15. The integral from 770 to infinity is then carried out 
analytically with the asymptotic formula (14). It turns out that an expansion in 
64 to 128 polynomials for each coordinate is sufficient to provide a numerically 
optimal approximation within the used precision. This illustrates the convenience 
of the disk coordinates in this context. The virial identity is then satisfied to the 
order of 10~^^. We plot the deviation of the sum of the integrals in (45) from zero 
for several values of A and 7 in Fig. 9. The drop in accuracy for strongly relativis- 
tic almost static situations (7 small and e close to 1) is again due to the almost 
degenerate Riemann surface. The lower accuracy in the case of strongly relativistic 
situations for 7 = 1 reflects the fact that the disk is shrinking to a point in this 
limit. To maintain the needed resolution one would have to use more polynomials 
in the evaluation of the virial-type identity which was not possible on the used 
computers. 
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Figure 9. The deviation from zero of the virial-type identity for 
7 = 0.7 and several values of e, and for e = 0.8 and several values 
of 7. 



5. Testing LOREME 

One purpose of exact solutions of the Einstein equations is to provide test-beds 
for numerical codes to check the quality of the numerical approximation. In the 
prc'vious sections we have established that the thcta-functional solutions can be 
numerically evaluated to the order of machine precision which implies they can be 
used in this respect. 

The code we are considering here is a C-|— I— library called LORENE [35] which was 

constructed to treat problems from rclativistic astrophysics such as rapidly rotating 
neutron stars. The main idea is to solve Poisson-type equations iteratively via 
spectral methods. To this end an equation as the Ernst equation (1) is written in 
the form 

(31) I^T = g{T,r,6,<i>), 

where spherical coordinates r, 6, (j) are used, and where Q is some possibly non- 
linear functional of and the coordinates. The system (31) is to be solved for T 
which can be a vector. In an iterative approach, the equation is rewritten as 

(32) /yFn+i=Q{J'n,r,e,cj>), n=l,2,.... 

Starting from some initial function JFqj in each step of the iteration a Poisson 
equation is solved for a known right-hand side. For the stationary axisymmctric 
Einstein equations which we are considering here, it was shown in [31] that this 
iteration will converge exponentially for small enough boundary data if the initial 
values are close to the solution of the equation in some Banach space norm. It 
turns out that one can always start the iteration with Minkowski data, but it is 
necessary to use a relaxation: instead of the solution Tn+i of (32), it is better to 
take a combination Tn+i = J^n+i + t^J^n with k €]0, 1[ (typically k ~ 0.5) as a new 
value in the source Qn+i to provide numerical stability. The iteration is in general 
stopped if ||jr„+i - J^„|| < 10-1°. 
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The Ernst equation (1) is already in the form (31), but it has the disadvantage that 
the equation is no longer strongly elliptic at the ergo-sphere where ^{S) = 0. In 
physical terms, this apparent singularity is just a coordinate singularity, and the 
theta-functional solutions are analytic there. The Ernst equation in the form (31) 
has a right-hand side of the form '0/0' for ^£ = which causes numerical problems 
especially in the iteration process since the zeros of the numerator and the denom- 
inator will only coincide for the exact solution. The disk solutions we are studying 
here have ergo-spheres in the shape of cusped toroids (see [13]). Therefore it is 
difficult to take care of the limit 0/0 by using adapted coordinates. Consequently 
the use of the Ernst picture is restricted to weakly relativistic situations without 
ergo-spheres in this framework. 

To be able to treat strongly relativistic situations, we use a different form of the 

stationary axisymmctric vacuum Einstein equations which is derived from the stan- 
dard 3 + 1-decomposition, see [16]. We introduce the functions v and via 

(33) e- = ^^, N,- ^"^^ 



where ae^^ is the gt^, component of the metric leading to the Ernst potential, 

see (37) in the appendix. Expressions for a in terms of thcta-functions are given 
in [13]. The vacuum Einstein equations for the functions (33) read 

(34) = ^p^e-^-^(Nl^ + Nl^), 

(35) AN^-^N^ = 4p(iV^,,(e2'^), + Ar^,f(e2-)c). 

By putting V = cos (f) we obtain the flat 3-dimensional Laplacian acting on V 
on the left-hand side, 

(36) AV = 4p{V,{e'np + Vde'nc)- 

Since the function e^'^ can only vanish at a horizon, it is globally non-zero in the 
examples we are considering here. Thus the system of equations (34) and (36) is 
strongly elliptic, even at an ergo-sphere. 

The disadvantage of this regular system is the non-linear dependence of the poten- 
tials 1/ and on the Ernst potential and a via (33). Thus we loose accuracy due 
to rounding errors of roughly an order of magnitude. Though we have shown in the 
previous sections that we can guarantee the numerical accuracy of the data for / 
and af to the order of lO"^**, the values for v and V are only reliable to the order 
of 10-13. 

To test the spectral methods implemented in LOREME, we provide boundary data for 

the disk solutions discussed above on a sphere around the disk. For these solutions 
it would have been more appropriate to prescribe data at the disk, but LORENE was 
developed to treat objects of spherical topology such as stars which suggests the 
use of spherical coordinates. It would be possible to include coordinates like the 
disk coordinates of the previous section in LORENE, but this is beyond the scope of 
this article. Instead we want to use the Poisson-Dirichlet routine which solves a 
Dirichlet boundary value problem for the Poisson equation for data prescribed at 
a sphere. We prescribe the data for v and on a sphere of radius R and solve 
the system (34) and (36) iteratively in the exterior of the sphere. If the iteration 
converges, we compare the numerical solution in the exterior of the sphere with the 
exact solution. 

Since spherical coordinates are not adapted to the disk geometry, a huge number of 
spherical harmonics would be necessary to approximate the potentials if R is close 
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to the disk radius. The limited memory on the used computers imposes an upper 
limit of 64 to 128 harmonics. We choose the radius R and the number of harmonics 
in a way that the Fourier coefficients in 9 drop below 10^^* to make sure; that the 
provided boundary data contain the related information to the order of machine 
precision. The exterior of the sphere where the boundary data are prescribed is 
divided in two domains, one from R to 2R. and one from 2R to infinity. In the 
second domain 1/r is used as a coordinate. For the (p dependence which is needed 
only for the operator in (36), 4 harmonics in are sufficient. 

Since LORENE is adapted to the solution of the Poisson equation, it is to be expected 

that it reproduces the exact solution best for nearly static situations, since the 
static solutions solve the Laplace equation. The most significant deviations from 
the exact solution are therefore expected for ^ = 0. For the case A = 3, we consider 
32 harmonics in on a sphere of radius R ~ 1.5. The iteration is stopped if 
||J>j+i — J>j < 5 * 10"^" which is the case in this example after 90 steps. The exact 
solution is reproduced to the order of 10~^^. The absolute value of the difference 
between the exact and the numerical solution on a sphere of radius 3 is plotted in 
Fig. 10 in dependence of 9. There is no significant dependence of the error on 9. 
The maximal deviation is typically found on or near the axis. As can be seen from 



maxdifftheta.eps 



Figure 10. Difference between the exact and the numerical solu- 
tion for A = 3 and 5 = for r = 3 in dependence on 6. 



Fig. 11 which gives the dependence on r on the axis, the error decreases almost 
linearly with 1/r except for some small oscillations near infinity. 

We have plotted the maximal difference between the numerical and the exact so- 
lution for a range of the physical parameters A and 6 in Fig. 12. As can be seen, 
the expectation is met that the deviation from the exact solution increases if the 
solution becomes more relativistic (larger e). As already mentioned, the solution 
can be considered as exactly reproduced if the deviation is below 10^^''. Increasing 
the value of 7 for fixed e leads to less significant effects though the solutions become 
less static with increasing 7. 
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Figure 1 1 . Difference between the exact and the numerical solu- 
tion for A = 3 and (5 = on the axis in dependence on r. 



gamma. eps 



Figure 12. Difference between the exact and the numerical solu- 
tion for 7 = 0.7 and several values of e, and for e = 0.8 and several 
values of 7. 



For (5 = 0, the ultra-relativistic limit A 4.629 . . . corresponds to a space-time with 
a singular axis which is not asymptotically flat, see [13]. Since LORENE expands all 
functions in a Galerkin basis with regular axis in an asymptotically flat setting, 
solutions close to this singular limit cannot be approximated. Convergence gets 
much slower and can only be achieved with considerable relaxation. For A = 4 and 
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(5 = we needed nearly 2000 iterations with a relaxation parameter oi k = 0.9. The 
approximation is rather crude (in the order of one percent). For higher values of A 
no convergence could be obtained. 

This is however due to the singular behavior of the solution in the ultra-relativistic 
limit. In all other cases, LORENE is able to reproduce the solution to the order of 
10~^^ and better, more static and less relativistic cases are reproduced with the 
provided accuracy. 



6. Conclusion 



In this article we have presented a scheme based on spectral methods to treat 
hyperelliptic theta-functions numerically. It was shown that an accuracy of the 
order of machine precision could be obtained with an efficient code. As shown, 
spectral methods are very convenient if analytic functions are approximated. Close 
to singularities such as the degeneration of the Riemann surface, analytic techniques 
must be used to end up with analytic integrands in the discussed example. 

The obtained numerical data were used to provide boundary values for the code 
LORENE which made possible a comparison of the numerical solution to the bound- 
ary value problem with the numerically evaluated theta-functions. For a large range 
of the physical parameters the numerical solution was of the same quality as the 
provided data. The main errors in LORENE are introduced by rounding errors in the 
iteration. This shows that spectral methods provide a reliable and efficient numer- 
ical treatment both for elliptic equations and for hyperelliptic Riemann surfaces. 
However, to maintain the global quality of the numerical approximation an analyti- 
cal understanding of the solutions is necessary in order to treat the non-analyticities 
of the solutions. 



Appendix A. Einstein equations and integral identities 

The Ernst equation has a geometric interpretation in terms of the stationary ax- 
isymmetric Einstein equations in vacuum. The metric can be written in this case 
in the Weyl-Lewis-Papapetrou form (see [27]) 

(37) ds'^ = Qabd-x^dx^ = -f{dt + ad(f>f + {e^^{dp'^ + dC^) + p^d(j)'^)/f, 

where p and ( are Weyl's canonical coordinates and dt and are the commuting 
asymptotically timelike respectively spacelike Killing vectors. 

In this case the vacuum field equations are equivalent to the Ernst equation (1) 
for the complex potential £. For a given Ernst potential, the metric (37) can be 
constructed as follows: the metric function / is equal to the real part of the Ernst 
potential. The functions a and k can be obtained via a line integration from the 
equations 



(38) = 2p 



and 

(39) '^ = ^^-^^wfw- 

This implies that a is the dual of the imaginary part of the Ernst potential. The 
equation (39) for k follows from the equations 

(40) R^0 = ^^{Sjp), a, /3= 1,2,3, 
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where R is the (three-dimensional) Ricci tensor corresponding to the spatial metric 
h = diag(e^'', e^'', p^). This reflects a general structure of the vacuum Einstein 
equations in the presence of a Killing vector. For the Ricci scalar one finds 

(41) -^e'''R=kpp + kcc. 

We denote by h the determinant of the metric h. 

The Komar integral [24, 34] of the twist of the timelike Killing vector = dt over 
the whole spacetime establishes the equivalence between the asymptotically defined 
ADM mass and the Komar mass rriK, 



(42) 2 / dV (r^b - IgabT^ ) =■ ^k, 



isk 



where the integration is carried out over the disk, where Ua is the normal at the 
disk, and where Tab is the energy momentum tensor of the disk given in [13]. In 
other words the ADM mass can be calculated either asymptotically or locally at 
the disk. 

To obtain an identity which does not involve only surface integrals, we consider as 
in [15] an integral over the trace of equation (40) for the Ricci-tensor, 



(43) R : 



2/2 • 



To avoid numerical problems at the set of zeros of /, the so-called ergo-sphere 

(sec [13] for the disk solutions studied here), we multiply both sides of equation (40) 
by f^. Integrating the resulting relation over the whole space- time, we find after 
partial integration 
(44) 

pl poo POO POO POO 

- / dppfkc+ / dp / da{pf)pkp+{pf)ckc) = dp dCpf{£p£p+£ck) 

Jo Jo J -oo Jo J -oo 

here the only contributions of a surface integral arise at the disk, since fc cx 1/r^ 
for r ^ oo and since the axis is regular {k vanishes on the axis). If we replace k 
via (39), we end up with an identity for the Ernst potential and its derivatives, 

/■I q l*oo /"OO 

- / dpp-'f{£p£c + £^£p) + - / / dpd(:p\£p{£l + £l) + £p{£l + £l)) 
Jo ^ Jo Jo 

pOC pOO 

(45) =2 / dpdCpf£c£c. 
Jo Jo 

This identity (as the identity given in [15]) can be seen as a generalization of the 

Newtonian virial theorem. The relation (45) coincides with the corresponding rela- 
tion of [15] only in the Newtonian limit. This reflects the fact that generalizations of 
a Newtonian result to a general relativistic setting are not unique. Our formulation 
is adapted to the Ernst picture and avoids problems at the crgo-sphercs, thus it 
seems optimal to test the numerics for Ernst potentials in terms of theta-functions. 
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